arXiv: 1507.01055v5 [cond-matsoft] 20 Sep 2016 


Antipolar ordering of topological defects in active liquid crystals 


Anand U. Oza^ and Jorn Dunkel^ 

^ Courant Institute of Mathematical Sciences, 251 Mercer Street, New York, NY 10012, USA 
^Department of Mathematics, Massachussetts Institute of Technology, 

77 Massachusetts Avenue, Cambridge, MA 02139-4307, USA 
(Dated: September 21, 2016) 

ATP-driven microtubule-kinesin bundles can self-assemble into two-dimensional active liquid crys¬ 
tals (ALCs) that exhibit a rich creation and annihilation dynamics of topological defects, reminiscent 
of particle-pair production processes in quantum systems. This recent discovery has sparked con¬ 
siderable interest but a quantitative theoretical description is still lacking. We present and validate 
a minimal continuum theory for this new class of active matter systems by generalizing the classical 
Landau-de Gennes free-energy to account for the experimentally observed spontaneous buckling of 
motor-driven extensile microtubule bundles. The resulting model agrees with recently published 
data and predicts a regime of antipolar order. Our analysis implies that ALCs are governed by the 
same generic ordering principles that determine the non-equilibrium dynamics of dense bacterial 
suspensions and elastic bilayer materials. Moreover, the theory manifests an energetic analogy with 
strongly interacting quantum gases. Generally, our results suggest that complex non-equilibrium 
pattern-formation phenomena might be predictable from a few fundamental symmetry-breaking and 
scale-selection principles. 


Active materials [T] assembled from intracellular com¬ 
ponents, such as DNA [2], microtubules and motor pro¬ 
teins I3H5] promise innovative biotechnological applica¬ 
tions, from microscale transport and medical devices [2] 
to artificial tissues [T] and programmable soft materi¬ 
als ms]. Beyond their practical value, these systems 
challenge theorists to generalize equilibrium statistical 
mechanics to far-from-equilibrium regimes mis]. Recent 
experimental advances in the self-assembly and manipu¬ 
lation of colloids with DNA-mediated interactions [271- 
[29] have stimulated theoretical analysis that may even¬ 
tually help clarify the physical principles underlying self¬ 
replication [30H32] and evolution in viruses [33U35] and 
other basic biological systems. Yet, despite some partial 
progress uniiiHiiiiiiMin], our conceptual understand¬ 
ing of active materials, and living matter in general, re¬ 
mains far from complete. We do not know whether, or 
under which conditions, ‘universality’ ideas [40] that have 
proved powerful in the description of equilibrium systems 
can be generalized to describe collective dynamics of ac¬ 
tive matter not just qualitatively but also quantitatively. 
This deficit may be ascribed to the fact that mathemat¬ 
ical models have been successfully tested against experi¬ 
ments in only a few instances laiiiisiiiTHia]. 

Recently discovered 2D active liquid crystal (ALC) 
analogs 0 EM] comprise an important class of non¬ 
equilibrium systems that allows further tests of general 
theoretical concepts m and specific models. ALCs are 
assemblies of rod-like particles that exhibit non-thermal 
collective excitations due to steady external jSlE] or in¬ 
ternal EE] energy input. At high concentrations, ALCs 
form an active nematic phase characterized by dynamic 
creation and annihilation of topological defects OlllIlT], 
reminiscent of spontaneous particle-pair production in 
quantum systems. This phenomenon was demonstrated 
recently um [48] for ATP-driven microtubule-kinesin 
bundles trapped in flat and curved interfaces. Moreover, 


these experiments [48] revealed an unexpected nematic 
ordering of topological defects which is unaccounted for 
in current theoretical models. Understanding the emer¬ 
gence of such topological super-structures is crucial for 
the development and control of new materials, as recently 
demonstrated for colloidal liquid crystals [494(51] . 

We here develop and test a closed continuum theory for 
dense ALCs by generalizing the higher-order scalar and 
vector theories of soft elastic materials [52] and bacte¬ 
rial fluids naini to matrix-valued fields. Specifically, 
we propose a modification of the commonly-adopted 
Landau-de Gennes (LdG) free energy to account for the 
inherently different microscopic buckling behaviors of 
passive and active LCs [5]. While bending is energeti¬ 
cally unfavorable in passive LCs and hence penalized by 
the LdG energy functional, kinesin-driven ALCs buckle 
spontaneously even at low concentrations due to the ex¬ 
tensile motor action (Fig. la). This experimental obser¬ 
vation [5] implies that the classical LdG framework is, by 
construction, ill-suited to describe experiments in which 
microtubule bundles are sheared relative to each other by 
motor proteins 0I171I18|. The inclusion of active bend¬ 
ing effects in the LdG functional yields a tensor version 
of the Swift-Hohenberg theory m of pattern formation. 
The resulting minimal model has only two dimensionless 
parameters, thus allowing a detailed comparison with re¬ 
cent experimental data [5] [48] . 

In addition to the traditional Q-tensor formulation, 
we present an equivalent complex scalar field represen¬ 
tation [HI El that manifests an analogy with a general¬ 
ized Gross-Pitaevskii theory ISSIEI of strongly coupled 
many-body quantum systems [574(60] . In the case of nor¬ 
mal dispersion, the celebrated LC-superconductor corre¬ 
spondence El ED has helped elucidate profound paral¬ 
lels between the smectic phase in passive LCs and the 
Abrikosov vortex lattices in type-II superconductors [62] 
[63] . The results below indicate that a similar analogy 
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may exist between ALCs and Bose-Einstein/Fermi con¬ 
densates with double-well dispersion [58H6Q] . suggesting 
that ALCs could offer insights into the dynamics of these 
quantum systems and vice versa. 


RESULTS 

Experimental conditions. Recent experiments [SJ 
[48] show that ATP-driven microtubule-kinesin bundles 
can self-assemble into a dense quasi- 2 D ALC layer at 
a surfactant-supported oil-water interface parallel to a 
planar solid boundary (Fig. [^). This ‘wet’ ALC was 
found to exhibit local nematic alignment of bundles, per¬ 
sistent annihilation and creation dynamics of topological 
defects [5], and remarkable nematic order of the defect 
orientations in thin layers [48| . Although a large number 
of unknown parameters has prevented detailed quanti¬ 
tative comparisons between theory and experiment, sev¬ 
eral recently proposed multi-order-parameter models of 
2 D ALC systems [HI |20l [Ml ES] were able to repro¬ 
duce qualitatively selected aspects of these observations, 
such as defect-pair creation and separation [65] . Despite 
providing some important insights, traditional models 
often do not account for three relevant details of the 
experiments (5] [48]. First, those models typically as¬ 
sume divergence-free 2 D ffuid ffow within the ALC layer, 
which is a valid approximation for isolated free-standing 
film experiments [ 66 ] but neglects ffuid exchange be¬ 
tween the 2 D interface and bulk in the ALC experi¬ 
ments (Fig. [^). Indeed, the surfactant-stabilized inter¬ 
face causes the microtubule-kinesin bundles to assemble 
into a quasi-2D layer, but places no such constraint on 
the ffuid. As is known for classical turbulence [67] [68] . 
small-scale energy input can trigger turbulent upward 
cascades in incompressible 2D ffow. Thus, topological 
defect dynamics in the current standard models may be 
dominated by artificially enhanced hydrodynamic mix¬ 
ing due to a simplifying 2D incompressibility assumption 
that is unlikely to hold under realistic experimental con¬ 
ditions [SIIIH]. Second, a relevant yet previously ignored 
effect is damping from the nearby boundaries, which 
may promote topological defect ordering. Third, as al¬ 
ready mentioned above, the commonly adopted standard 
LdG free-energy functional does not account for motor- 
driven spontaneous buckling [5] of microtubule bundles 
(Fig . [^) , which is one of the key differences between pas¬ 
sive and active LCs (Z. Dogic, private communication). 
To overcome such limitations and achieve a quantitative 
description of the experiments la EH] , we next construct a 
closed continuum theory for ALCs described by a nematic 
tensor field The theory accounts for the different 

buckling behaviors of passive and active LCs and builds 
on a self-consistent hydrodynamic closure condition. 

Theory. Traditional multi-field models [64l ES] aim to 
describe the 2 D nematic phase of a dense ALC suspension 
by coupling the filament concentration c(t, r) and the ne¬ 
matic order tensor Q(t^r) to an incompressible 2 D flow 


field v{t^r) that satisfies V • v = 0 in the interface plane 
r = (x^y). The nematic order parameter 5'(t,r) is pro¬ 
portional to the larger eigenvalue of Q, and the filaments 
are oriented along the corresponding eigenvector, or di¬ 
rector n(t,r). To construct an alternative closed-form 
theory for the symmetric traceless 2 x 2 -tensor field Q, 
we start from the generic transport law 

dtQ+ V ■ (vQ) - k[Q,uj] =-^ (1) 

where uj = [\/v — {\/v)^]/2 is the vorticity tensor, 
[A^B] = AB — BA the commutator of two matrices and 
^[Q] = JcfvF an effective free energy. Focussing on 
dense suspensions as realized in the experiments [5] [48] , 
we neglect fluctuations in the microtubule concentra¬ 
tion, Vc = 0. A derivation of the advection term 
V • {vQ) from the probability conservation laws underly¬ 
ing generic advection-diffusion models is outlined in the 
Supplementary Information. It is important, however, 
that V • {vQ) ^ V 'VQ when V-v 7 ^ 0, which is typically 
the case when fluid can enter and leave the interface. 
Combining LdC theory [69] with Swift-Hohenberg the- 
ory [ 53 ], we consider the effective non-equilibrium free- 
energy density (Supplementary Information) 

F = IV - |(VQ)2 + ^(vvg)2| (2) 

with a, 6 > 0 for the nematic phase. Assuming 72 can 
have either sign, ultraviolet stability requires 74 > 0. For 
72 < 0 , F penalizes bending and buckling, as is appro¬ 
priate for passive LCs and possibly ‘dry’ shaken nemat¬ 
ics [111 [ 45 ] , forcing the system dynamics towards a ho¬ 
mogeneous nematic ground-state manifold. By contrast, 
motor-induced spontaneous buckling [5] (Fig. [^) of 
kinesin-driven ALCs demands 72 > 0, and co nseque ntly 
patterns of characteristic wavelength A ^ A/ 74/72 be¬ 
come energetically favorable, as shown below. 

The requirement 72 > 0 for ALCs has an intrinsi¬ 
cally microscopic origin, as the ALC assembly consists of 
microtubules that grow against each other and sponta¬ 
neously buckle due to the motor-induced extensile shear 
dynamics of adjacent bundles (Fig. [^). To explain this 
important point in more detail, let us recall that passive 
LCs are modeled using 72 < 0 , as the corresponding term 
in the free energy penalizes variations in Q and thus in¬ 
hibits spatial inhomogeneities at damping rate 72 < 0 
in Fourier space. Microscopically, the alignment dynam¬ 
ics of two rod-like passive LC molecules roughly corre¬ 
sponds to a time-reversal of the ALC pair-interaction 
image sequence in Fig. implying that a correspond¬ 
ing ALC systems would develop buckling instabilities at 
growth rate 72 > 0 . Additional empirical support for 
this theoretical picture comes from a comparison of the 
experimentally observed length scales in ALC systems: 
The microtubule-kinesin bundles realized in the ALC ex¬ 
periments [5] [48] are approximately 10-30 /im in length 
(Fig. lb in [5]), which is on the same scale as both the 
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FIG. 1: (a) Image sequence showing spontaneous buckling of a microtubule bundle (dashed) caused by extensile ATP-driven 
motor activity, adapted with permission from Fig. Ic in Ref. [5]. Time interval 11.5s, scale bar 15 /im. For a passive or 
contractile bundle, one would expect straightening instead of bending, approximately corresponding to a time-reversal of the 
depicted sequence, (b) Schematic of the experimental setup reported in Ref. O [48], not drawn to scale. A thin oil him 
(thickness ~ 3/xm) separates a 2D ALC layer (~ 0.2 — 1.0/xm) at the oil-water interface from a solid glass cover. Liquid can 
be exchanged between the ALC layer and bulk huid, resulting in compressible 2D interfacial how that is strongly damped by 
the nearby no-slip glass boundary and the viscous oil layer, (c) Extensile 2D dipole how in the interface as predicted by the 
overdamped closure condition 0 for D > 0 and Q = (A, 0; 0, — A) with A = exp(—r^). The central horizontal bar indicates the 
unit director axis, and background colors the nematic order parameter S ^ X. 


spontaneous buckling wavelength (Fig. Ic, d in [5|, and 
Fig.[^) and the typical separation distance between de¬ 
fects (bottom panels of Fig. 3d in [5], reproduced in 
Fig. below). That is, there exists no substantial scale 
separation between the buckling microscopic constituents 
and emergent ALC dynamics. 

Similar buckling phenomena are generically observed 
in many systems that are subjected to external or inter¬ 
nal stresses, for example in elastic hlms and sheets [5^(70] 
and in geometrically conhned cellular networks [111172]. 
The ALCs experience an ehFective compressive stress due 
to the extensile ’growth’ of the hlament pairs in a con¬ 
hned geometry, which arises from their motor-induced 
shear dynamics. Application of such a compressive stress 
leads to buckling of the network’s constituents [niizi]. 
It has been shown that out-of-plane buckling of an elas¬ 
tic sheet due to an effective compressive stress may be 
quantitatively modeled by a Swift-Hohenberg-type equa¬ 
tion with 72 > 0 in the corresponding free energy [52]. 
We expect the same to be true for the in-plane buckling 
of microtubules conhned to a planar interface, and thus 
analyze here an ehFective theory that incorporates this 
spontaneous motor-induced buckling phenomenologically 
through 72 > 0. 

Hydrodynamic closure. To obtain a closed Q- 
model, we relate the 2D how held v to Q through the 
linearly damped Stokes equation [731174] 

— r]\/^v uv = —CV • Q (3) 

where r] is the viscosity and the rhs. represents active 
stresses [I0l|6l] with C > 0 for extensile ALCs (Fig.[^). 
A pressure term does not appear in Eq. because 
the interfacial how is not assumed to be incompress¬ 
ible and concentration huctnations are neglected. The 
h'-teim in the force balance (§ has been used to model 
interfacial damping in other contexts, such as surfactant 


membranes on a solid substrate m, and accounts for 
friction from the nearby no-slip boundary in the Hele- 
Shaw m approximation (Fig. [^). In the overdamped 
regime /rj ^ 1, we deduce from Eq. (0 the closure 
condition 

i; = -DV • Q , D = C^/v. (4) 

Equation 0 is conceptually similar to closure con¬ 
ditions proposed previously for active polar hlms 
m- Importantly, Eq. 0 predicts divergent interfacial 
how, V • V 7 ^ 0, and hence huid transport perpendicular 
to the interface wherever VV : Q 7 ^ 0. Inserting 0 
into Q yields a closed Q-theory in which periodic di¬ 
rector patterns corresponding to local minima of the free 
energy can become mixed by self-generated interfacial 
how. 

Complex representation ALC-quantum anal¬ 
ogy. The traditional characterization of 2D nematic 
order in terms of the symmetric traceless 2 x 2 matrix 
held Q = (A,/x;/x,—A) is redundant, for only two real 
scalar helds A(t,r) and iJi{t^r) are needed to specify 
the nematic state at each position r = {x^y). To 
obtain an irreducible representation [M] [54] we de- 
hne the complex position coordinate z = x veloc¬ 
ity held v{t,z) =u^iw and complex order parameter 
'ip{t^z) = A + i/x, such that S = 2|'0|. In terms of the 
Wirtinger gradient operator d=^{dx—idy)^ the 2D 
Laplacian takes the form = Add and the closure con¬ 
dition Q reduces to x; = —2Dd^l^. Denoting the real and 
imaginary parts of an operator O by and ^{O}, 

Eqs. 0 and 0 may be equivalently expressed as 

dt'ip + A^^ = -^ (5) 

where the self-advection operator is given by 

= -4D5R{(57) + (5V’)5} + 4/tF>i9{aV} (6) 
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and the free energy Q['ip^'ip] = J dz G has the density 

G = + 72'0(459)'0 + 74'0(495)^'0. (7) 

For 72 < 0 and 74 ^ 0, Eq. 0 reduces to the energy den¬ 
sity of the Gross-Pitaevskii mean-field model [55l[56] for 
weakly interacting boson gases. Historically, this limit 
case has been crucial [sa EH for elucidating the anal¬ 
ogy between the smectic phase of passive LCs and the 
Abrikosov flux lattice in type-II superconductors [62l [63] . 
For 72,74 > 0, Eq. 0 effectively describes double-well 
dispersion EU, as recently realized for quasi-momenta in 
spin-orbit-coupled Bose-Einstein condensates isHiisni and 
Fermi gases [59] . This fact establishes an interesting con¬ 
nection between dense ALCs and strongly coupled quan¬ 
tum systems: when self-advection is negligible ( 1 ) —^ 0 ), 
the fixed point configurations of Eq. 0 coincide with the 
‘eigenstates’ of generalized Gross-Pitaevskii models that 
incorporate wavelength selection. 

Stability analysis. The qualitative model dynam¬ 
ics is not significantly altered for moderate values of k 
(Movies SI and S 2 ), so we neglect the commutator term 
by setting k = 0 from now on (see Supplementary In¬ 
formation for hi > 0). To understand the properties of 
Eqs. ([^-(|^ when self-advection is relevant, we perform 
a fixed point analysis of the rescaled dimensionless equa¬ 
tion (Supplementary Information) 

dt'ip — 4:D -h {d^)d} = 

Q - V’ - 72(4aa)tA - i4.ddf^p, ( 8 ) 

by focussing on the uniform state which 

corresponds to a nematic order parameter value S = 1 
and homogeneous director angle 0 relative to the x-axis. 
Considering wave-like perturbations pj = ^ e(t)e*^'^ 

with |e| <C 1 and extensile ALCs with D > 0, 
one finds that ip^ is unstable when 72 > 0 (Supple¬ 
mentary Information). For subcritical self-advection, 
D < Dc = 2(—72 + y /72 + 2 ), the dominant instability is 
driven by modes with wavenumber \k\ = y^ 72 /( 274 ), 
suggesting the formation o f stripe p atterns with typ¬ 
ical wavelength A « ^ 87 r^ 74 / 72 . By contrast, 

for supercritical advection, D > Dc^ the most un¬ 
stable mode propagates perpendicular to the direc- 

tor, (fc*,(/)*) = (y(272 +D)/(474),6' + (7r/2)), suggest- 
ing the possibility of transverse mixing. 

Phase diagram. To investigate the nonlinear dynam¬ 
ics of Eq. 0 , we implemented a Fourier pseudospectral 
algorithm with modified Runge-Kutta time-stepping Hg 
(Methods) and so evolved the real and imaginary parts 
of ' 0 (t, z) in time for periodic boundary conditions in 
space. A numerically obtained ( 72 ,D)-phase diagram 
for random initial conditions confirms the existence of 
a turbulent nematic phase if active self-advection is suf¬ 
ficiently strong (Fig. |^,b; Movie SI). Ordered configu¬ 
rations prevail at low activity (Fig. ^,c,d; Movies S3, 
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FIG. 2: Phase diagram obtained from simulations of Eq. Q 
for one particular set of random initial conditions showing 
the emergence of turbulent nematic states for supercritical 
active self-advection. (a) We observe convergence to defect- 
free stripes (blue, panel c, Movie S3), long-lived static and 
oscillatory defect lattice solutions (green, panel d, Movie S4), 
oscillatory defect creation and annihilation events (black. 
Movies S6 and S7), and chaotic dynamics (red, panel b, 
Movie SI). The white line indicates the analytical estimate 
Dc = 2 (—72 + ^7! + 2) for the transitions between ordered 
and chaotic states, (b-d) Examples of the states identified 
in a with -defects (black) and +|-defects (white). Panel 
d highlights the antipolar ordering of +|-defect orientations 
(red bars); see also Fig. and Movie S5 for a realization of 
this state in a ~ lOx larger simulation domain. 


S4 and S 6 ). Although the ground-states of the free 
energy 0 are in general not homogeneous, the criti¬ 
cal curve separating the two regimes is i n fair agree¬ 
ment with the estimate Dc = 2(—72 + y ^72 + 2 ) from 
linear stability of the homogenous state (white line in 
Fig. [^). For subcritical values of the advection parame¬ 
ter D^ we observe either defect-free ground-states or long- 
lived lattice-like states exhibiting ordered defect config¬ 
urations (Fig. 0 ,d). Regarding the subsequent compar¬ 
ison between theory and experiment, it is important to 
note that the lattices are also found in simulations with 
a large domain (Fig. Movie S5). These spatially peri¬ 
odic states generally exhibit antoolar long-range order¬ 
ing of +|-defects (FigjM; Fig.^) accompanied by vor¬ 
tex flow lattices (Fig. [3^. Numerical free-energy calcula¬ 
tions show that defect-free states (Fig.|^) typically have 
slightly lower energies than the lattice states (Fig. 0 ), 
leaving open the possibility of a very slow decay of the 
latter. However, regardless of whether such lattice states 
are extremely long-lived metastable or truly stable states, 
these simulation results confirm that antipolar ordering 
of +|-defect pairs can persist over experimentally rele- 
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FIG. 3: Vortex lattice states with antipolar long-range ordering of nematic defects, see also Movie S5. (a) Long-lived nematic 
order parameter field with periodically aligned — |-defects (black) and +|-defects (red), (b) Line integral convolution (LIC) 
plot of the corresponding director field as a proxy for microtubule-bundle patterns, (c) LIC plot of the corresponding fluid 
velocity field, color-coded by normalized vorticity, demonstrates the formation of a vortex flow lattice. Dimensionless simulation 
parameters are D = 0.25, 72 = 1.875. 



FIG. 4: Defect-pair creation and propagation in experiment and theory, (a) Experimentally observed dynamics of a defect pair, 
spontaneously produced by buckling and subsequent fracture of filaments; adapted with permission from Fig. 3d in Ref. [5]. 
Scale bar 20 /xm, time lapse 15 s. (b) Line integral convolution (LIC) plot of the director fields showing the spontaneous creation 
and propagation of a defect-pair in a simulation of Eq. ^ for D = 1.5, 72 = 1. As in the experiments, +|-defects (yellow) 
generally move faster than — |-defects (light blue), cf. Eig. 


vant time-scales. 

Theory vs. experiment. To test our theory sys¬ 
tematically against existing experimental data [3 EH] , we 
analyze defect-pair dynamics, global defect ordering and 
defect statistics in the turbulent nematic phase. 

Spontaneous defect-pair creation and subsequent prop¬ 
agation, as reported in recent ALC experiments [5] and 
observed in our simulations, are compared in Fig. In 
the experimental system [ 5 ], a (+^, -defect pair is 
created when fracture along incipient crack regions [20] 
becomes energetically more favorable than buckling. Af¬ 
ter creation, the +|-defect moves away rapidly whereas 
the position of the —^-defect remains approximately 
fixed for up to several seconds (Fig. Ef)- We note that 
an asymmetry in the speeds of topological defects has 
also been observed in passive liquid crystals [77H79] . Our 
simulations of the minimal model defined in Eq. ® ac¬ 
curately reproduce the details of the experimentally ob¬ 
served dynamics (Fig. Ef; Movies SI and S2). 

Another striking and unexplained experimental obser¬ 
vation [48] is the emergence of orientational order of +|- 
defects in thin ALC layers (Fig. [^. Using the setup 
illustrated in Fig. [^, recent experiments [48] demon¬ 
strated nematic alignment of +|-defects in thin ALC 
layers of thickness h ^ 250 nm (Fig. [^), whereas thicker 


ALC layers with h ^ 1 /xm showed no substantial ori¬ 
entational order on large scales (Fig. |^). To inves¬ 
tigate whether our theory can account for these phe¬ 
nomena, we tracked defect positions (Methods) and 
defect orientations di = V • (5(^i)/|V • (5(^i)| [HO] in 
simulations for different values of the advection param¬ 
eter D = C/z^, since Brinkman-type scaling arguments 
suggest that D increases with the ALC layer thickness, 
D (X Iju (X with p G [1,2]. For weakly super¬ 
critical advection, D > Dc^ we find that Eq. ^ pre¬ 
dicts robust antipolar alignment of +l-defects (Fig. ^). 
In our simulations, this ordering decreases as the effec¬ 
tive mixing strength D increases (Fig. |^), consistent 
with the experimental results [48] for thicker ALC layers 
(Fig. [^). Similar ordering was observed in simulations 
that incorporated alignment of the nematic field near 
the horizontal boundaries of the simulation box (Sup¬ 
plementary Information; Fig. S5). To quantify the de¬ 
gree of orientational order, we recorded the distances dij 
between all +1-defect pairs (lj) as well as their rela¬ 
tive orientation angles Oij = cos“^(d^ • dj) G [0,7r]. The 
resulting pair-orientation distributions p{0\r)^ and polar 
and nematic correlation functions, P(r) = {di • dj)r and 
N{r) = 2{{di'djY)r — l, are shown in Fig.j^ with (•)^ de¬ 
noting an average over pairs of defects separated by a dis- 
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FIG. 5: Strong and weak antipolar ordering of +|-defects as 
(a,b) observed in experiments [48] and (c,d) predicted by our 
theory based on 2D simulations with periodic boundary condi¬ 
tions. Light-blue markers: — |-defects. Yellow markers: +|- 
defects. Red bars: orientation of +|-defects. (a) +|-defects 
in thin ALC films (thickness h ~ 250 nm) show strong nematic 
alignment, (b) +|-defects in thicker ALC films (h ~ 1 fim) 
are more disordered. (c,d) For weak effective hydrodynamic 
coupling D, simulations show antipolar ordering, which is in¬ 
hibited for larger values of D. The average number of defects 
in the full simulation box is approximately (c) 240 and (d) 
350. Figures a and b kindly provided by S. DeCamp and Z. 
Dogic. 


tance r. The local maxima in the orientation distribution 
at 6> = 0 and 0 = ir signal antipolar ordering (Fig. [^), 
which is also reflected in the oscillatory behavior of the 
polar and nematic correlation functions (Fig.[^,d). The 
diminished intensity of the local maxima for larger val¬ 
ues of D indicates that enhanced hydrodynamic mixing 
reduces orientational order (Fig. |^,d). 

Lastly, we test the theoretically predicted defect statis¬ 
tics against a separate experimental data set kindly pro¬ 
vided by DeCamp et al. (private communication). Since 
our simulations are performed in dimensionless units, 
there is freedom to choose a characteristic lengthscale 
/o and timescale to- To relate theory and experiments, 
we determine (/o, ^o) such that the joint mean speed and 
mean lifetime of ±^-defects match the experimental val¬ 
ues V = 6.6 /im/s and r = 52.8 s. After fixing these global 
scales, we can compare details of the speed and lifetime 
distributions (Fig. [^. To this end, we first locate the 
‘best-fit’ simulation parameters in the ( 72 , 1 ^)-parameter 
space explored in the phase diagram (Fig. [^). This pro¬ 
cedure identifies 72 = = 1.75 as the best-match 

parameters, although nearby parameter values and sim¬ 
ulations with n = 1 produce fits of similar quality, cor¬ 
roborating the robustness of the model (Fig. S4). For 


— ^-defects, we find adequate agreement between exper¬ 
iment and theory for speed and lifetime probability den¬ 
sity functions (PDFs), as evident from Fig. [^,b. For 
+ ^-defects, simulation results also agree well with the 
experimental measurements (Fig. [^,d), but one notices 
two systematic differences. First, while the peak heights 
of the PDFs agree within a few percent, experimen¬ 
tally measured speed values for +|-defects are on av¬ 
erage slightly larger than theoretically predicted values 
(Fig. ^). Second, simulation data predict a miniscule 
tail-fraction of long-living -defects not detected in the 
experiment (Fig. IT)- In addition, based on the experi¬ 
mental density estimate of 30 defects/mm^ [48], we find 
that the defect density at any given time in the ‘best-fit’ 
simulation is ^ 2 . 3 x lower than in the experiments. As 
discussed below, such deviations can be explained plau¬ 
sibly by specific model assumptions. Taken together, the 
above results confirm that the minimal model defined by 
Eq. provides a satisfactory qualitative and quantita¬ 
tive description of the main experimental results (5] |48] . 





+1/2-Pair distance r/ro +1/2-Pair distance r/ro 

FIG. 6: Increasing activity and film thickness decreases 

antipolar ordering in simulations. (a,b) Maxima of the nu¬ 
merically obtained local pair orientation PDFs p{0ij\r) signal 
antipolar local ordering of +|-defects as they are separated 
by the typical defect-lattice spacing. The defect distance r 
is specified in units of the mean nearest-neighbor distance 
ro between +|-defects. (c,d) Polar P{r) and nematic N{r) 
correlation functions for D = 1.5 (red) and D — 3 (blue). 
Increasing the effective hydrodynamic coupling D leads to 
stronger mixing and hence decreases nematic order, which is 
corroborated by the nematic correlation length being ~ 40% 
shorter for D — ?> than for D = 1.5 (panel d). This is also 
reflected by the diminished intensity of the maxima in panel 
b relative to panel a. The simulation parameters correspond 
to those given in Fig. [^,d. 
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FIG. 7: Quantitative comparison of defect statistics between 
predictions of Eq. ^ and experimental data |38], using the 
parameter estimation procedure described in the text. For 
— I-defects, both (a) speed distribution and (b) lifetime dis¬ 
tribution agree well, (c) For +|-defects, experimentally mea¬ 
sured speed values are slightly larger, as our model assumes 
a strongly overdamped limit, (d) Simulations with periodic 
boundary conditions (Movie S 8 ) predict a low-probability tail 
of large lifetimes which is not visible in the experiment, likely 
due to its restricted field of view or additional noise. Di¬ 
mensionless simulation parameters D = 1.75 and 72 = 1 
translate into the following dimensional values: a = 0.08 s“^, 
b = 0.32s“^, D = 1791/im^/s, 72 = 1024/im^/s, 74 = 
3.28 X lO^/xm^/s. The numbers n± reflect the number of ±|- 
defects tracked, and the simulation domain contained ~ 130 
defects at any given time. 


DISCUSSION 

Pattern-formation mechanism. Equations 
and 0 epitomize the idea of ‘universality’ in spatio- 
temporal pattern formation, as known from Swift- 
Hohenberg-type scalar field theories [53l l 8 T] . The free- 
energy expressions contain the leading-order terms of 
generic series expansions in both order-parameter space 
and Fourier space, consistent with spatial and nematic 
symmetries. When considering passive systems with a 
preference for homogenization (72 < 0 ), it usually suf¬ 
fices to keep only the quadratic gradient terms. By 
contrast, for pattern forming systems, the coefficient in 
front of the lowest-order gradient contribution can change 
sign [ 52 I [ 53 ] , and one must include higher-order deriva¬ 
tives to ensure stability. We hypothesize that the sign 
change of 72 is directly related to the motor-induced 
buckling of microtubule bundles (Fig. [^), an effect that 
is not captured by the standard LdG free-energy for 
passive liquid crystals. In a few select cases, expres¬ 
sions of the form 0 and 0 can be systematically de¬ 


rived [52l |53l [82]. Generally, one can regard the free- 
energy expansion 0 as an effective field theory whose 
phenomenological parameters can be determined from 
experiments. This approach has proved successful for 
dense bacterial suspensions ns E] and now also for 
ALCs, indeed suggesting some universality in the for¬ 
mation and dynamics of topological defects in active sys¬ 
tems. 

Nematic defect order. Although ( 72 ,!^) are varied 
as independent effective parameters in the simulations, 
they are likely coupled through underlying physical and 
chemical parameters. For example, it is plausible that a 
change in ATP-concentration or film-thickness would af¬ 
fect both 72 and D. The parameter D can also be inter¬ 
preted as an effective Reynolds number. In our numer¬ 
ical exploration of the ( 72 , D)-parameter space, we ob¬ 
serve for subcritical advection D either long-lived lattice¬ 
like states exhibiting nematically aligned -defects or 
defect-free ground-states (Figs. Movies S3, S4, S5, 
S 6 and S7). Ordered defect configurations correspond to 
local minima or saddles in the free-energy landscape and 
have only slightly higher energy than defect-free states 
(Fig. §). When the activity ( is sufficiently large that 
advection is marginally supercritical, D > Dc^ chaotic 
system trajectories spend a considerable time in the 
vicinity of these metastable lattice states, which provides 
a physical basis for the orientational order of defects in 
thin ALC films [48]. For D ^ Dc, the ALC system can 
access a wider range of high-energy states, leading to in¬ 
creased disorder in the defect dynamics. Although the 
strongly turbulent regime D ^ Dc requires high time- 
resolution and is thus difficult to realize in long-time sim¬ 
ulations, the inhibition of nematic defect order at larger 
values of D is evident from the reduced peak heights in 
Fig. and d. 

Defect statistics. The systematic speed deficit in 
Fig. likely reflects the overdamped closure condi¬ 
tion H, which suppresses the propagation of hydrody¬ 
namic excitations. Since flow in a newly created defect 
pair generally points from the — ^ to the -defect, the 
minimal model ^ can be expected to underestimate the 
speeds of + ^-defects. This effect could be explored in 
future experiments through a controlled variation of the 
thickness and viscosity of the oil film (Fig.[^). The low- 
probability tail of long-living +|-defects in the simula¬ 
tion data (Fig. 1^) may be due to the fact that they can 
be tracked indefinitely in the simulations but are likely 
to leave the finite field of view in the experiments. In 
the future, the minimal theory presented here should be 
extended systematically by adding physically permissible 
extra terms [83] to the free-energy, explicitly simulating 
the full hydrodynamics in Eq. 0 , or incorporating addi¬ 
tional terms into Eq. 0 that account for the interaction 
between the nematic field and the induced flow {65] [83] . 

Future extensions. The minimal model formulated 
in Eqs. 0 and 0 can be systematically extended to 
improve further the quantitative agreement between ex¬ 
periment and theory. For instance, one may append to 

























the right-hand side of Q an additional fourth-order lin¬ 
ear term of the form (VV)+ [V • (V • Q)] [84], 0+ de¬ 
noting the symmetric traceless part of the operator O. 
Such a term only affects the high-wavenumber damp¬ 
ing at order and thus is not expected to alter sig¬ 
nificantly the results obtained here. We also note that 
extra terms coupling the nematic field to the induced flow 
may be added to Eq. p), an example being SE, where 
E = (1/2) \Vv (Vv) ' is the symmetrized strain rate 
tensor [251165]. Our above analysis neglected such sec¬ 
ondary hydrodynamic effects in the interest of construct¬ 
ing a minimal mathematical theory capable of captur¬ 
ing key experimental observations. Moreover, this sim¬ 
plification may be justified on the physical basis that 
steric interactions and motor-induced buckling are ex¬ 
pected to dominate over flow-alignment effects at high 
microtubule densities^. The effect of microtubule bend¬ 
ing may be enhanced by appending a hydrodynamic in¬ 
teraction term proportional to E+ [25] , since the closure 
condition v = —DV • Q implies that E^ = —{D/2)AQ, 
which augments the bending term —j 2 ^Q in Eq. 
However, the scale separation between the experimen¬ 
tally observed filament buckling wavelength and the flow 
structures in the isotropic phase at low microtubule con¬ 
centrations (see Eig. Id in Ref. [5]) suggests that such 
hydrodynamic effects play a secondary role. A natural 
next step would be to derive systematically the bending 
term from a microscopic model of motors and filaments 
as introduced in Ref. [20]. Einally, it will be worthwhile 
to attempt constructing a fully 3D theory for the ALCs 
and fluid and subsequently project on the quasi-2D inter¬ 
face, although additional assumptions are required then 
to obtain a closed 2D system of equations. 


CONCLUSIONS 

Recent experimental and theoretical studies showed 
that fourth-order PDE models for scalar and vector 
fields provide an accurate quantitative description of 
surface-pattern formation in soft elastic materials [52] 
and orientational order in dense bacterial fluids [13 SI]. 
Here, we have generalized these ideas to matrix-valued 
fields describing soft active nematics. The above anal¬ 
ysis demonstrates that a generic fourth-order Q-tensor 
model can shed light on experimental observations in 2D 
ALCs [48], including the emergence of orientational or¬ 
der of topological defects. Physically, the higher-order 
generalization Q becomes necessary because the com¬ 
monly adopted LdG free-energy, which was designed to 
describe passive liquid crystals, does not account for the 
experimentally observed spontaneous buckling of motor- 
driven ALCs [5]- More generally, the fact that three 


vastly different soft matter systems can be treated quan¬ 
titatively in terms of structurally similar higher-order 
PDEs [T3|4TJ[52] promises a unified mathematical frame¬ 
work for the description of pattern formation processes 
in a broad class of complex materials. In addition, the 
free-energy analogy [54] between dense ALCs and gen¬ 
eralized Gross-Pitaevskii models suggests that the self¬ 
organization principles m of mesoscopic active matter 
and microscopic quantum systems [FTHGQ] could be more 
similar than previously thought. 
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METHODS 

Numerical solver defect tracking. To sim¬ 
ulate Eq. we implemented a numerical algorithm 
that evolves the real and imaginary parts A(t,r) and 
/i(t,r) of the complex order parameter 7 / in time for 
periodic boundary conditions in space. The algorithm 
solves Eq. pseudospectrally in space using N£ = 256 
or Ni = 512 lattice points in each direction and a simu¬ 
lation box of size L = Gtt (Eig. 2 and 4) or L = ISir (Eig. 
3, 5-7). Spectral analysis shows that N£ = 256 is gener¬ 
ally sufficient to resolve the fine-structure of the numer¬ 
ical solutions (Supplementary Eig. S6). The algorithm 
steps forward in time using a modified exponential time- 
differencing fourth-order Runge-Kutta method [76] with 
time step At < Simulations were initialized with 

either a single defect pair or random field configurations 
{A(0, r), /i(0, r)}. Defects are located at the intersections 
of the zero-contours of A and /i, their positions tracked by 
implementing James Munkres’ variant of the Hungarian 
assignment algorithm [86] (Supplementary Information). 


^ In active nematics of low to intermediate density, concentration 
fluctuations can trigger additional instabilities m- 
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